FinalPaper html version

Introduction and research question

Life-Expectancy has constantly increased over the last centuries, and it has never been so high as nowadays. Medicine has taken great strides forward, eradicating and finding cures for a lot of mortal diseases. Having said that, there are still other big challenges that we have to face in medicine and one of those is for sure cancer. Cancer is a leading cause of death all over the world, accounting for almost 10 milion deaths in 2020. This is the reason why we thought that would be very interesting to analyze the principal factors that, according to state-of-the-art scientific research, cause this disease. We also wanted to divide these factors in some macro-groups of variables: environmental, life-style related and socio-economical/demographic.

Our analysis (purely statistical) aims to verify if it is indeed possible using data about the selected risk factors to predict the share of people with cancer in a given year and in a given country. Moreover, we tried to improve our all-inclusive model by finding out the best model out of the above mentioned factors

Dataset, Data Cleaning and Data Visualization

Source

Our research started by searching out on some reliable websites what are the leading causes of cancer and dividing them into the macro-categories. To do that we adopted the WHO official website and we selected the following risk factors:

• ENVIRONMENTAL: Air pollution
• LIFE-STYLE RELATED: Alcohol consumption, Tobacco consumption, Obesity
• SOCIO-ECONOMICAL/DEMOGRAPHIC: GDP, Share of old people

(IMPORTANT DISCLAMER: of course the factors defined under “Socio-Economical/Demographic” cannot be considered to cause cancer, but as we will see they turn out to be very useful to model cancer all over the world and over time)

Description of the variables

We built our dataset selecting by merging different csv files all taken from the same source, that is called “Our world in data”. This is by itself collecting datasets from verified sources and the precise reference to each factor is given in the sitography at the end of the paper. Here we provide a complete description of the variable we decided to use in our model.

Cancer: Share of total population with any form of cancer, measured as the age-standardized percentage. This share has been age-standardized assuming a constant age structure to compare prevalence between countries and through time.

Air_pollution: Population-weighted average level of exposure to concentrations of suspended particles measuring less than 2.5 microns in diameter (PM2.5). Exposure is measured in micro grams of PM2.5 per cubic meter (µg/m³).

Alcohol: Average per capta consumption of alcoholic beverages, measured in kilograms per year. Data is based on per capta food supply at the consumer level, but does not account for food waste at the consumer level

GDP per capta: Measured in constant international-$

Obesity: Obesity is defined as having a body-mass index (BMI) equal to or greater than 30. BMI is a person’s weight in kilograms divided by his or her height in metres squared

Old age Dependency: This is the ratio of the number of people older than 64 relative to the number of people in the working-age (15-64 years). Data are shown as the proportion of dependents per 100 working-age population.

Smoking: Share of population who smoke every day

Constructing the table

We started by loading all the different csv files. Then we took the column we were interested in and we glued them together. After some other little manipulation (see the R code for that) here is how part of the final table looks like (the complete table with the N.A.s results in a dimension of 6468x11):

##We start by loading some useful libraries that we will need during this process

library(readr) #Used to read the csv
library(tidyverse) #Contains a lot of useful packages to clean the datas
library(gganimate) #To create animated plots
library(ggthemes) #To select some nice themes 
library(zoo) #To fill the N.A. values 
library(grid) #To display plots into a grid
library(gridExtra) #To display plots into a grid
library(ggpubr)#To display plots into a grid

##We now load the row csv files what we need to merge 

Cancer <- read_csv("Datasets/Cancer.csv")
Air_pollution <- read_csv("Datasets/Air pollution.csv")
Alcool <- read_csv("Datasets/Alcool.csv")
GDP_per_capta <- read_csv("Datasets/GDP per capta.csv")
Obesity <- read_csv("Datasets/Obesity.csv")
Old_age_dependency_ratio <- read_csv("Datasets/Old age dependency ratio.csv")
Smoking <- read_csv("Datasets/Smoking.csv")
Population <- read_csv("Datasets/population-since-1800.csv")
countryContinent <- read_csv("Datasets/countryContinent.csv")

#We proceed by taking the explanatory variable data set (Cancer) and adding the covariates as columns other columns, to do that we use the use the libraby "dplyr" as it makes it very easy to do and readable

Continent = countryContinent %>% select(country, continent)
Continentt = rename(Continent, Entity = country)
Cancer1 = left_join(Cancer, Air_pollution, by = c("Entity", "Code", "Year"))
Cancer2 = left_join(Cancer1, Alcool, by = c("Entity", "Code", "Year"))
Cancer3 = left_join(Cancer2, GDP_per_capta, by = c("Entity", "Code", "Year"))
Cancer4 = left_join(Cancer3, Obesity, by = c("Entity", "Code", "Year"))
Cancer5 = left_join(Cancer4, Old_age_dependency_ratio, by = c("Entity", "Code", "Year"))
Cancer6 = left_join(Cancer5, Smoking, by = c("Entity", "Code", "Year"))
Cancer7 = left_join(Cancer6, Continentt, by = "Entity")
Cancer_final = Cancer7

rm("Cancer1", "Cancer2", "Cancer3", "Cancer4", "Cancer5", "Cancer6", "Cancer7")
rm("Air_pollution", "Alcool", "Cancer", "GDP_per_capta", "Obesity", "Old_age_dependency_ratio", "Smoking", "Continent", "Continentt", "countryContinent")

#We now proceed by renaming the columns. This will make all the manipulations of the dataframe more readable:

colnames(Cancer_final)[4] <- "Cancer"
colnames(Cancer_final)[5] <- "AirPoll"
colnames(Cancer_final)[6] <- "Alcool"
colnames(Cancer_final)[7] <- "GDP"
colnames(Cancer_final)[8] <- "Obesity"
colnames(Cancer_final)[9] <- "Age"
colnames(Cancer_final)[10] <- "Smoking"

Cancer_final

write.csv(Cancer_final, "Cancer_final.csv", row.names=FALSE, quote=FALSE) #to save the file into our environment

Data Visualization

Now we will proceed with some plotting. Again, for the detailed process for making this graphs please refer to the R script. You may also want to check our our file in Html format, in which we made some nice animation to be able to display the evolution in time of the variables.

We start by displaying the variable we are interested to explain, that is the share of population with any form of cancer. We chose to plot the most recent year in our data frame, that is 2017. As we can clearly see form the plot below the countries in which are suffering more (in percentage) of cancer are USA, Canada, Australia and Greenland. These states are followed by European countries. Asian and African countries are the ones with relatively less cases of cancer in percentage. Here we could argue that a big downside of this picture is the fact that different countries have different ability in the process of screening. As the WHO writes “Late-stage presentation and lack of access to diagnosis and treatment are common, particularly in low- and middle-income countries”. We will further discuss this limitation in the dedicated section. (Countries for which data were not available are omitted from the map, but those are luckily few)

#We read the file needed
share_of_population_with_cancer <- read_csv("Cancer.csv")  #read the file 
share_of_population_with_cancer = data.frame(lapply(share_of_population_with_cancer, function(x) {
  gsub("United States", "USA", x) #correction needed for the mapping 
}))
colnames(share_of_population_with_cancer)[1] <- "region" #Give the same name for the country in the two datasets
colnames(share_of_population_with_cancer)[4] <- "cancer" #Dropping the complicated name given by the data frame

# Create the actual map that we will use
#We will use the ggplo2 package for that

mapdata = map_data("world") #Used to take the latitude and longitude of the countries 
mapp <- mapdata %>%
  left_join(share_of_population_with_cancer, by = "region") #merge the two datasets
mapp2 = mapp %>% 
  filter(!is.na(mapp$cancer)) #dealing with NA
mapp2$cancer <- as.numeric(mapp2$cancer) #changing the value to numeric for plotting
mapp2$Year <- as.numeric(mapp2$Year) #changing the value to numeric for animation

# Building the map
#Here we create a static map
ggplot( data = mapp2, aes(long, lat, group = subregion)) +
  geom_map(
    aes(map_id = region),
    map = mapdata,
    color = "black", fill = "gray30", size = 0.3
  ) +
  geom_polygon(aes(group = group, fill = cancer), color = "black") +
  scale_fill_gradient2(low = "blue", high = "red" , midpoint = 1.0355583) +
  theme_minimal() +
  transition_time(Year)
  labs(
    title = "Share of population with cancer in {round(frame_time}",
    x = "", y = "", fill = ""
    ) +
  theme_bw()
  shadow_wake(wake_length = 0.05)
animate(graphPoll, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Share with cancer world.gif")

-> This is not the file you would get running the code, but its static version, My computer just doesn’t have enough computational power to produce this animation and shuts down when I try to create the animation : (

Now we are going to explore the dataset by plotting each selected factor against the Cancer one to start having an idea of some possible relation. Because we selected them as “risk factors” we would expect to see some positive correlation going on. To be consistent with what we have done above we again plot the graphs relative to the year 2017. We also used some libraries to make plots nicer by coloring the points according to the continent and by changing the size according to the population of the country. The result is the following:

#Here we correct the N.A.s by filling the closest non empty value (this is needed for plotting)
#and we will not do any analysis with these fake data
Cancergraph = left_join(Cancer_final, Population, by = c("Entity", "Code", "Year"))
Cancergraph$Smoking = na.locf(Cancergraph$Smoking)
Cancergraph$AirPoll = na.locf(Cancergraph$AirPoll)
Cancergraph$Alcool = na.locf(Cancergraph$Alcool)
Cancergraph$GDP = na.locf(Cancergraph$GDP, fromLast = TRUE)
Cancergraph$Age = na.locf(Cancergraph$Age)
Cancergraph$Obesity = na.locf(Cancergraph$Obesity)
Cancergraph$`Population (historical estimates)` = na.locf(Cancergraph$`Population (historical estimates)`)

#Here we create and save each plot, This chunk of code will not be evaluated every time because is pretty heavy and time/CPU consuming. The results (trust me) are included at the end
graphPoll = Cancergraph %>%
  ggplot(aes(x = AirPoll, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs Air Pollution", 
       x = "Air pollution exposure ", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphPoll, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Airpoll Visual Animation.gif")


graphAlco = Cancergraph %>%
  ggplot(aes(x = Alcool, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs Alcool consumption", 
       x = "Alcool consumed", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphAlco, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Alcool Visual Animation.gif")


graphGDP = Cancergraph %>%
  ggplot(aes(x = GDP, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs GDP", 
       x = "GDP per capta", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphGDP, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("GDP Visual Animation.gif")


graphObes = Cancergraph %>%
  ggplot(aes(x = Obesity, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs Obesity", 
       x = "Share of Obese people", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphObes, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Obesity Visual Animation.gif")


graphAge = Cancergraph %>%
  ggplot(aes(x = Age, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs Age", 
       x = "Old people ratio", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphAge, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Age Visual Animation.gif")


graphSmoke = Cancergraph %>%
  ggplot(aes(x = Smoking, y = Cancer, color = continent, size = `Population (historical estimates)`)) +
  geom_point(alpha = 0.9, stroke = 0) +
  scale_size(range = c(2,12), guide = "none") +
  scale_color_brewer(palette = "Set2") +
  labs(title = "Cancer incidence vs Smoking", 
       x = "Number of Smokers", 
       y = "Share of people with Cancer",
       color = "Continent",
       caption = "Source: Our World in Data") +
  transition_time(Year) +
  labs (subtitle = "Year:{round(frame_time)}") 
 shadow_wake(wake_length = 0.05)
animate(graphSmoke, height = 500, width = 800, fps = 30, duration = 15, end_pause = 60, res = 100)
anim_save("Smoking Visual Animation 3.gif")

#And eventually here it is the animations:
knitr::include_graphics("Airpoll Visual Animation.gif")

knitr::include_graphics("Alcool Visual Animation.gif")

knitr::include_graphics("GDP Visual Animation.gif")

knitr::include_graphics("Obesity Visual Animation.gif")

knitr::include_graphics("Age Visual Animation.gif")

knitr::include_graphics("Smoking Visual Animation 3.gif")

Our hypothesis is confirmed by these graphs: evidently there’s some sort of positive correlation, but still not strong as we imagined beforehand. We can see that age is the factor which presents a stronger correlation visually.The only problem is the one related to Air pollution. Our explanation for that is that will be provided in the section “Limitations of the study”.

We also have some outliers, which are Canada and US (The two orange points that on each graph are above the others). Being cancer causes such a complicated and intricate thing, is actually quite difficult to really understand why those two countries have such an high number of cases. We can suppose that these two countries perform a more intense screening activity, but this goes beyond the scope of this study. We will just delete these two countries from our data set. Having said that, the result that follow will be not significantly changed by this action.

Multivariate Linear Regression

We are now going to perform a multivariate linear regression, trying to examine the relationship between some explanatory variable and one response variable. As mentioned before, we are going to take into consideration as explanatory variables all the factors selected to help us trying to explain part of the variance of the share of people affected by cancer worldwide. The table used is considering approximately 150 countries for approximately 20 years. In this way we will be able to do not miss some trend or to do not being deceived from some outliers years/countries or random trends.

Now, let’s proceed with performing the actual multiple linear regression. After that, we want to see if the assumptions required to state the statistical significance are met: we’ll check for the normality and for the homoscedasticity of the residuals.

library(olsrr)

# Delete the N.A.s
mydata = read.csv("Cancer_final.csv")
data5 = na.omit(mydata)

# Delete USA and Canada from the data set
data3 = subset(data5, Entity != "USA" || Entity != "Canada")

#Normality of the residuals
model = lm( Cancer~ Smoking + Obesity + GDP + Age + AirPoll + Alcool, data = data3)

#Perform Testing for normality of residuals
ols_test_normality(model)

## -----------------------------------------------
##        Test             Statistic       pvalue  
## -----------------------------------------------
## Shapiro-Wilk              0.722          0.0000 
## Kolmogorov-Smirnov        0.1671         0.0000 
## Cramer-von Mises         216.6397        0.0000 
## Anderson-Darling         46.5541         0.0000 
## -----------------------------------------------

#plot histogram of residuals
kop = residuals.lm(model)
hist(kop, breaks = 50)

#Check graphically using the function plot
plot(model)

Let’s analyze each of the test we performed above and see which conclusion we can take from what we observed. For what concerns normality of the residuals, we can see that the p-value of the tests we performed is below the 0.05 significance level for all of them, thus we can reject the null hypothesis that the errors are normally distributed. This is exactly as we imagined: the data set is made by more or less 1000 data points, thus violating this assumption was something that we expected for. This is due to the fact that with this many data points it’s easy for the testing methods to find evidence against normality. Trying to see the situation graphically, the situations seems a bit better. The qq-plot shows that the data fit not that badly the diagonal line, exceptions made for a little deviation near the top-right. Furthermore, analyzing another graph we see that plotting the residuals against the fitted values approximates more or less a straight line. The analysis of the residuals is enough satisfying, but if we want to be strict we will say that we are not able to take out any significant conclusions from the experiment we will carry on.

We now perform and interpret the result of the multiple linear regression itself.

summary(model)

## 
## Call:
## lm(formula = Cancer ~ Smoking + Obesity + GDP + Age + AirPoll + 
##     Alcool, data = data3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71319 -0.14060 -0.02478  0.08123  2.64302 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.020e-01  4.013e-02   5.035 5.74e-07 ***
## Smoking     -2.872e-03  1.371e-03  -2.095 0.036418 *  
## Obesity      2.647e-03  1.414e-03   1.871 0.061588 .  
## GDP          1.666e-05  9.009e-07  18.497  < 2e-16 ***
## Age          5.045e-02  2.250e-03  22.419  < 2e-16 ***
## AirPoll     -4.296e-03  6.737e-04  -6.377 2.83e-10 ***
## Alcool       1.150e-03  3.259e-04   3.528 0.000439 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3053 on 941 degrees of freedom
## Multiple R-squared:  0.8243, Adjusted R-squared:  0.8232 
## F-statistic:   736 on 6 and 941 DF,  p-value: < 2.2e-16

From just these results, we can see various things: the value of the Adjusted R-Squared is very satisfying, and then this model fits pretty well the data, telling us there is a quite good relationship with the covariates we indicated: our factors explain a 82% variance in the share of population with cancer! We are actually considering the Adjusted R-Squared since we are performing a multivariate linear regression, and thus the number of factors we include in our model is relevant and we want to penalize bigger models. In fact we know that the simply value of R-Squared can just increase when adding more factors, even if those are completely random ones. In addition, all the variables seem to explain pretty well the share of cancer cases:almost all of the p-values are statistically significant! The only exception is given by the “Obesity” explanatory variable.

Model Selection (Testing)

Now that we discussed the linear model and the underlying characteristics, we want to do some basic model selection to see if there are other models (possibly simpler) that are able to explain the dependent variable well enough. To do that we will leverage on two different approach to model selection. The first is the one provided by testing and the other the one provided by information theory. We will compare the results and see whether they coincide or not.

We start by importing the data we will need for the model selection and loading some useful libraries. We will then create our linear model using all the covariates and the dataset without the missing values (we will use the same model obtained above with the regression)

library(olsrr) #This stands for Ordinary Least SquaRe Regression and it has the tools needed
model = lm( Cancer ~ Smoking + Obesity + GDP + Age + AirPoll + Alcool, data = data3)
summary(model)

## 
## Call:
## lm(formula = Cancer ~ Smoking + Obesity + GDP + Age + AirPoll + 
##     Alcool, data = data3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71319 -0.14060 -0.02478  0.08123  2.64302 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.020e-01  4.013e-02   5.035 5.74e-07 ***
## Smoking     -2.872e-03  1.371e-03  -2.095 0.036418 *  
## Obesity      2.647e-03  1.414e-03   1.871 0.061588 .  
## GDP          1.666e-05  9.009e-07  18.497  < 2e-16 ***
## Age          5.045e-02  2.250e-03  22.419  < 2e-16 ***
## AirPoll     -4.296e-03  6.737e-04  -6.377 2.83e-10 ***
## Alcool       1.150e-03  3.259e-04   3.528 0.000439 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3053 on 941 degrees of freedom
## Multiple R-squared:  0.8243, Adjusted R-squared:  0.8232 
## F-statistic:   736 on 6 and 941 DF,  p-value: < 2.2e-16

We begin by trying three different testing methods, called respectively “Step-up/Forward”, “Step-down/Backward” and “Step-both”. The fist two act in a dual way with respect to the second and the last is in a way able to overcome the downsides of the first two models. Without entering too much into the details this greedy methodology starts with no(all) predictors in the model and iteratively add(subtract) the one that is more statistically significant in contributing to the model. The procedure automatically stops when no other variables can be added to gain statistical significance. In practice what we do at each step is to perform a t-test for each candidate variable and see which one has the lowest p-value (if there is one). In this setting we kept a significance level of 0.05 both for entering and for exiting the model.

The Step-up method yields:

frwfit <- ols_step_forward_p(model, pent = 0.05)
frwfit

## 
##                             Selection Summary                              
## --------------------------------------------------------------------------
##         Variable                  Adj.                                        
## Step    Entered     R-Square    R-Square      C(p)        AIC        RMSE     
## --------------------------------------------------------------------------
##    1    Age           0.7156      0.7153    579.6198    896.8510    0.3875    
##    2    GDP           0.8108      0.8104     71.6307    512.4933    0.3162    
##    3    AirPoll       0.8210      0.8204     19.0720    462.0425    0.3078    
##    4    Alcool        0.8231      0.8224      9.5664    452.6010    0.3061    
## --------------------------------------------------------------------------

The selection led us to add to the model the “Age”, “GDP”, “Airpoll” and “Alcool”. We now try the symmetrical method,

The Step-down method yields:

bkwfit <- ols_step_backward_p(model, prem = 0.05)
bkwfit

## 
## 
##                           Elimination Summary                            
## ------------------------------------------------------------------------
##         Variable                  Adj.                                      
## Step    Removed     R-Square    R-Square     C(p)       AIC        RMSE     
## ------------------------------------------------------------------------
##    1    Obesity       0.8237      0.8227    8.5024    451.5307    0.3057    
##    2    Smoking       0.8231      0.8224    9.5664    452.6010    0.3061    
## ------------------------------------------------------------------------

As you can see we got the same result here: we rejected the “Obesity” and the “Smoking” cofactor. Notice that even if the methods seem completely specular, it’s not granted that they’ll lead us to the same outcome.

To do one last cross-check using testing methods we use a mixed method to allow variables that were added in previous steps to be excluded in succeeding steps (this wasn’t allowed in the simple forward/step-up method)

The “Step-both” method yields:

bothfit <- ols_step_both_p(model, pent = 0.05, prem = 0.05)
bothfit

## 
##                               Stepwise Selection Summary                               
## --------------------------------------------------------------------------------------
##                      Added/                   Adj.                                        
## Step    Variable    Removed     R-Square    R-Square      C(p)        AIC        RMSE     
## --------------------------------------------------------------------------------------
##    1      Age       addition       0.716       0.715    579.6200    896.8510    0.3875    
##    2      GDP       addition       0.811       0.810     71.6310    512.4933    0.3162    
##    3    AirPoll     addition       0.821       0.820     19.0720    462.0425    0.3078    
##    4     Alcool     addition       0.823       0.822      9.5660    452.6010    0.3061    
## --------------------------------------------------------------------------------------

But this again yields the same result as the one we got with the step forward method. As well as before, the fact that we obtained the same result with all the three approaches is a very likable. This seems to suggest a model of four variables that is cutting “Smoking” and “Obesity”, given that these two risk factor are removed by all three of our greedy methods.

We now try all together different methods based on Information criterion. These work in such a way that they penalize big models that otherwise would be selected by methods such as LSE and MLE. The math underneath this approach is absolutely non trivial and explaining it goes beyond the scope of this research. To implement this we will use the function “best_subset” That is returning the best model containing from one to six covariates. We therefore limit our comperison to the comparison between all the best subsets models that we can have.

allsub2 <- ols_step_best_subset(model)
allsub2

This seems to suggest a different results. All of these criteria are suggesting that the best model is the one with all six the explanatory variable with the exception of the SBC method, that again is suggesting the model with four variables as before. This discrepancy can be explained by the fact that, at least for the Akaike Information Criteria that we studied in class, the underlying constant C could dominate. This method is in fact more suitable for models with different dimensions |d|. (Notice also the the difference in AIC is very low: \(\Delta AIC = 2.59\) ) With small models AIC can easily prefer the complete model as the gain in “explained variance” obtained by adding a variable is more than the penalty introduced by that extra parameter.

We conclude this section by saying that the best model after this analysis is the one with four explanatory variable, obtained with the three step-wise method

Limitations of the study and Conclusions

Based on the results of our study, we can conclude that our last model explains quite well the variance in the data and that we collected about cancer and that some factors (for sure including the ones selected by us) are more influential than others in explaining the share of people suffering of cancer in a given country. Relatively older and richer countries (the two things often overlap) present many more cancer cases in percentage than the others and the same we can conclude of countries with a consumption of alcohol and cigarettes above the mean. However, our study has a lot of limitations that we have to indicate. - First, as we had to construct the data set by our own, we encountered some difficulties: we joined some data together into a table, but obviously we had some missing data. It’s not easy even just to collect the data in the same way for all the countries for more than 25 year for all those factors! - Secondly, we have to consider the difference in data availability that is present among the countries. To give an example: it’s much easier to find a data set about a given topic regarding the USA with respect to one regarding Burkina Faso. Besides that, we have taken into considerations just six factors to try to explain cancer incidence, but as we know it is a much more difficult research question: genetical heritage, other factors of concerning lifestyle, other environmental elements, access to the healthcare system, prevention and treatment offered by the country and many more things play a role in developing a cancer or not. - Moreover, it is not an ‘immediate’ disease: factors of risk sum up over the years and so do progress in medicine. This ends up in a complicated and intricate situation to analyze. Diving into some details, we also had a quite surprising result: the correlation between cancer cases and air pollution seems strongly negative, as the plot at the beginning shows. It is difficult to explain such a result, but it may be that since the worst quality of air is on average in less developed countries(eg. Nigeria, Pakistan, India…), it may be that other factors as the demographic distribution of the population or GDP, that typically obtain low values in those countries , have a stronger influence on the actual number of cancer cases, ending in this counter-intuitive result. The division that we made at the beginning seems to suggest that the “non-risk factors” are the one better suited to describe the distribution of cases worldwide, but again that’s what we wanted for this research: a merely descriptive approach. We conclude by saying that another interesting improvement to this research could be to trying to repeat the same analysis but just focusing on everyday habits and maybe by including the interaction effects within the possible covariates that we will choose in order to provide tips for an healthier routine.

->Hope many of you got the StatQuest tribute<-

SITOGRAPHY AND REFERENCES

The sources used for the datas are:

• Cancer -> Global Burden of Disease Collaborative Network. Global Burden of Disease Study 2017 (GBD 2017) Results. Seattle, United States: Institute for Health Metrics and Evaluation (IHME), 2018

• Air_pollution -> World Development Indicators - World Bank

• Alcool -> Food and Agriculture Organization of the United Nations (FAO) (2017)

• GDP_per_capta -> World Development Indicators - World Bank

• Obesity -> World Health Organization (WHO)

• Old_age_dependency -> World Development Indicators - World Bank

• Smoking -> Nationally representative sources, survey data